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The gamma-ray population of millisecond pulsars (MSPs) detected by the Fermi Large Area Telescope 
(LAT) has been steadily increasing. A number of the more recent detections, including PSR J0034— 0534, 
PSR J1939-I-2134 (B1937+21; the first MSP ever discovered), PSR J1959-I-2048 (B1957-I-20; the first black 
widow system), and PSR J2214-(-3000, exhibit an unusual phenomenon: nearly phase-aligned radio and gamma- 
ray light curves (LCs). To account for the phase alignment, we explore geometric models where both the radio 
and gamma-ray emission originate either in the outer magnetosphere near the light cylinder (-Rlc) or near 
the polar caps (PCs). We obtain reasonable fits for the first three of these MSPs in the context of "altitude- 
limited" outer gap (alOG) and two-pole caustic (alTPC) geometries. The outer magnetosphere phase-aligned 
models differ from the standard outer gap (OG) / two-pole caustic (TPC) models in two respects: first, the 
radio emission originates in caustics at relatively high altitudes compared to the usual low-altitude conal radio 
beams; second, we allow the maximum altitude of the gamma-ray emission region as well as both the minimum 
and maximum altitudes of the radio emission region to vary within a limited range. Alternatively, there also 
exist phase-aligned LC solutions for emission originating near the stellar surface in a slot gap (SG) scenario 
("low-altitude slot gap" (laSG) models). We find best-fit LCs using a Markov chain Monte Carlo (MCMC) max- 
imum likelihood approach |3C| . Our fits imply that the phase-aligned LCs are likely of caustic origin, produced 
in the outer magnetosphere, and that the radio emission may come from close to Rlc- We lastly constrain 
the emission altitudes with typical uncertainties of ~ 0.3/Jlc- Our results describe a third gamma-ray MSP 
subclass, in addition to the two (with non-aligned LCs) previously found [50l| : those with LCs fit by standard 
OG / TPC models, and those with LCs fit by pair-starved polar cap (PSPC) models. 



1. INTRODUCTION 

The first pulsar catalog released by Fermi Large 
Area Telescope (LAT) included 46 gamma-ray pulsars 
0, 8 of which were millisecond pulsars (MSPs) 
Currently, there are > 20 gamma-ray MSPs [1^ and 
> 70 gamma-ray pulsars in total [4J|. The discov- 
ery of PSR J0034-0534 ^ revealed it to be the 
first MSP to have (nearly) phase-aligned radio and 
gamma-ray light curves (LCs). This rare phenomenon 
has only been observed for the Crab pulsar (32j . 
However, this behavior has now also been observed 
for PSR J1939-H2134 (B1937-1-21), PSR J19594-2048 
(B1957+20) [24i], and PSR J2214-I-3000 41], and more 
MSPs will be added to this subclass. 

1.1. Traditional Emission Models 

Two classes of pulsar models have been used to de- 
scribe high-energy (HE ) pulsar emission. In polar cap 
(PC) models [ij, [l^, primary electrons are ejected 
from the neutron star (NS) surface and accelerated 
along curved magnetic field lines, producing curva- 
ture radiation gamma rays. Thermal X-rays may also 
be upscattered to gamma-ray energies. Subsequently, 



these gamma rays are converted into electron-positron 
pairs via magnetic pair production in the intense mag- 
netic fields close to the stellar surface (at radius i?Ns)- 
In addition, a slot gap (SG) 0, [s^l may form along the 
last open magnetic field lines of the pulsar magneto- 
sphere in the absence of pair creation along those lines. 
This corresponds to a two-pole caustic (TPC) geom- 
etry [T6j which may extend from the stellar surface 
up to near the light c ylin der (at radius Rlc)- Outer 
gap (OG) models [lol.l43| represent the second model 
class. In these models, HE radiation is produced along 
the last open field lines above the null charge sur- 
face (NCS) where the Goldrcich- Julian charge density 
changes sign. The narrow gaps in both the OG and 
TPC models require screening of the electric field par- 
allel to the local magnetic field, and therefore presup- 
poses copious pair production. Lastly, HE LCs were 
also modeled in the context of OG and TPC models 
in a force-free magnetic field geometry, proposing a 
separatrix layer model close to i?LC • 

1.2. Formation of Caustics 

HE photons escaping from the magnetosphere are 
subject to two relativistic effects: their traveling direc- 
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tion is aberrated due to the large corotation velocity, 
while their arrival time at the observer is determined 
by their emission height, due to the finite speed of 
light. Lastly, it is assumed that these photons are 
emitted tangent to the local magnetic field lines in 
the co-rotating frame. The combination of these three 
effects result in the formation of caustics, i.e., the ac- 
cumulation of photons in narrow phase bands (36j . 
These caustics manifest themselves as bright peaks in 
the observed pulse profiles. 

1.3. MSP Models 

Due to their much lower surface dipole magnetic 
field strengths, MSPs have been thought to have 
pair-starved magnetospheres, where the magnetic pair 
multiplicity is not high enough to screen the accel- 
erating electric field in the open volume above the 
PC [iK,!!^. In this case, primary electrons are acceler- 
ated up to very high altitudes above the full PC, while 
pair formation is suppressed. This model is called 
the pair-starved polar cap (PSPC) model [s^, H^, an 
extension of the traditional PC model. MSP LCs 
and spectra have been modeled using this frame- 
work 111 m Ei]. Alternatively, MSP spectra and 
energetics have also been modeled in the context of 
an OG model [s^, [s^]- An annular gap model [l5| can 
furthermore reproduce the main characteristics of the 
gamma-ray LCs of three MSPs, although this model 
does not attempt to model the nonzero phase offsets 
between the gamma-ray and radio profiles. 

1.4. MSP Subclasses 

The first 8 i^erTni-detected gamma-ray MSPs have 
been modeled jS^I • Two distinct MSP subclasses were 
found: those whose LCs are well fit by a standard 
OG or TPC model, and those whose LCs are well 
fit by a PSPC model (with these fits being mutually 
exclusive). These models yielded the correct radio-to- 
gamma phase lags when the radio emission was mod- 
eled as a cone beam at lower altitude. Such fits im- 
plied that MSPs have screened magnetospheres with 
large amounts of pairs available, as these conditions 
are needed to set up the gap structure presupposed by 
the OG / TPC models. Small distortions of the dipole 
magnetic field causing offsets of the PC may provide 
a mechanism for enhancing pair creation, even in low- 
spin-down pulsars ^26]. This paper discusses a third 
sublcass of MSPs: those having phase-aligned radio 
and gamma-ray LCs. 

1.5. Motivation for Caustic Radio 
Emission 

In contrast to the first two MSP sublcasses, the 
near phase- alignment of the gamma-ray and radio LCs 



of MSPs in the third subclass argues for overlapping 
emission regions. These co-located emission regions 
may occur at high altitudes, so that the radio emis- 
sion will be subject to the same relativistic effects as 
the gamma-ray emission, as described in Section 11.21 
Discoveries of new gamma-ray MSPs exhibiting phase- 
aligned LCs therefore motivate the investigation of 
high-altitude [111 caustic radio emission. 

A second argument motivating caustic radio emis- 
sion comes from investigating the beaming properties 
of normal pulsars and MSPs detectable using blind 
searches on gamma-ray data as well as radio data [il] . 
The relative number of gamma-ray to radio pulsars 
for each of these 'gamma-ray-selected' and 'radio- 
selected' samples implies that radio and gamma-ray 
beams must have comparable sky coverage of ~ 47r sr 
for pulsars with high spin-down luminosities (Sj-ot), 
but radio beams should shrink for pulsars having lower 
values of i^rot- The radio emission for high-Sj-ot pul- 
sars should therefore originate in wide beams at a sig- 
nificant fraction of -Rlc- One should however bear 
in mind that LCs resulting from radio and gamma- 
ray caustics would generally be nearly phase-aligned 
(although small phase differences could result if the 
radio and gamma-ray emission regions are at different 
altitudes). Caustic radio emission is therefore plausi- 
ble for young pulsars with nearly aligned LCs. Radio 
caustics may however be more common in the case of 
the MSPs, as there are many more examples of MSPs 
with phase-aligned LCs. 

1.6. Modeling Phase-Aligned LCs 

We investigate the possibility of reproducing phase- 
aligned radio and gamma-ray LCs using "altitude- 
hmited" OG / TPC models (alOG / alTPC) in which 
we limit the extent of the emission regions (Sec- 
tion 13. ip vs. a low-altitude SG (laSG) model (Sec- 
tion By modeling the LCs of PSR J0034-0534, 
PSR J1939-h2134, and PSR J1959-I-2048, we can infer 
values for the magnetic inclination and observer angles 
a and C (Section [5]), and also constrain the emission 
altitudes. 



2. BACKGROUND ON SELECTED MSPs 
WITH PHASE-ALIGNED LCs 

2.1. PSR J0034-0534 

PSR J0034— 0534 was discovered using the Parkes 
radio telescope 6]. It follows a circular orbit around 
a low-mass companion {Hubble Space Telescope ob- 
servations revealed an optical white dwarf companion 
with a mass of about Q.2Mq; [8]). PSR J0034-0534's 
period of P = 1.87 ms implies a rotational age of Tc = 
P/2P - 10 Gyr, dipolar surface field of Ba ~ 10^ G, 
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and spin-down power of i^iot ^ 2 x 10^'* erg Q, 
typical among the radio MSP population. It is also 
relatively close, lying at 0.5 kpc [1^. No rotating 
vector model (RVM) fits exist for PSR J0034-0534 
pGj because there is no detected linear polarization. 
A < 3(7 detection of an X-ray source 0.2" from the 
pulsar position by XMM-Newton was reported [ssj ]. 
while EGRET obtained a 3a flux upper limit above 
100 MeV [IS exceeding the recent Fermi flux mea- 
surements [3|| by an order of magnitude. 



2.2. PSR J1 939+21 34 (B1 937+21) 

PSR J1939-h2134 is the first MSP ever discov- 
ered 0. It has a period P — 1.558 ms, and re- 
mains one of the fastest-spinning MSPs discovered. 
This MSP has a very high spin-down luminosity of « 

-4x10^ G, 
and lies at 
RXTE oh- 



Outer Gap 



\{yi& g-^.g g-i^ surface magnetic field of B\ 
and characteristic age of T^ 
a distance of cZ = 7.7 ± 3. 
servations 




ll| revealed a double-peaked X-ray profile 



with phase separation of about half a rotation, cl osely 
aligned with the phases of the giant radio pulses 3l| , 
but slightly lagging the radio peaks. The EGRET 
3(7 upper limit to the unpulsed flux above 100 MeV 
was 15.1 X 10~^ cm s~^ ^8j]. Pulsations with a sig- 
nificance well above 5ct have now been detected from 
PSR J1939+2134 by Fermi LAT [H. 



2.3. PSR J1 959+2048 (B1 957+20) 
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Figure 1: Schematic diagram of the emission layers in the 
OG (panel a) and TPC (panel b) geometries. The 
magnetic axis is indicated by pi, and the spin axis by f1. 
The two concentric circles indicate the limiting minimum 
and maximum emission radii imposed for the 
altitude-limited models. 



whether MSPs produce enough pairs to facilitate a 
force-free magnetosphere. 

We furthermore use a photon emission rate that is 
constant along magnetic field lines in the corotating 
frame, and treat relativistic effects (i.e., aberration 
of photon directions and time-of-flight delays; Sec- 
tion [L2]) consistently to first order in r/i?LC (with 
r the radial distance from the NS center). We lastly 
include the Lorentz transformation (a second-order ef- 
fect in r/i?Lc) of the local magnetic field between lab 
and corotating frames 0- 



PSR J1959-l-2048wasthe first "black widow" pulsar 
discovered. It is in a nearly circular eclipsing binary 
orbit, ablating its low-mass tidally-locked companion 
star [23, El]. Lying at a distance of d ^ 2.5 kpc, its 
period P = 1.607 ms and intrinsic P is 8x 10~^^ [11] 
lead to a spin-down luminosity of 7 x 10'^'' erg s~^, 
surface magnetic field of i?o ~ 10^ G, and character- 
istic age of Tc ^ 3 Gyr. XMM-Newton observations 
revealed a phase dependence of the X-ray emission 
on the binary orbital period [29j . although no pul- 
sations were detected at the pulsar spin period P. A 
~ 4(7 pulsed X-ray signal have now been observed from 
PSR J1959-I-2048, with the X-ray peaks seemingly in 
close alignment with the radio peaks [23|. Significant 
pulsed gamma-ray emission has also been detected by 
Fermi LAT [H]. 



3. GEOMETRIC PULSAR MODELS 

As in in our previous work |50l] , we assumed a 
retarded vacuum dipole magnetic field as the basic 
structure of the pulsar magnetosphere [H, [13] ■ In 
the case of young pulsars, this field may actually be 
closer to the force-free solution [i^ , but it is not clear 



3.1. High-altitude Gamma-ray and Radio 
Emission: alOG and alTPC Models 

We use the same framework as previously [50], 
but the radio emission region is extended in alti- 
tude. We free the minimum and maximum radii of 
the radio (i?Jni„ and i?Jnax) ^^d the maximum radius 
of the gamma-ray (i?2aax) emission regions, and re- 
strict the emission gaps' extent to a cylindrical radius 
Pmax < 0.95i?Lc- Importantly, we do not use an axis- 
centered radio conal model, but investigate radio pho- 
tons coming from an OG / TPG-like structure. In the 
alOG radio models, the minimum radius is actually 
max i?NCs}j so it is a function of magnetic az- 

imuth (p and co-latitude 9 when i?NCS > -Rmin- Here, 
-Rncs(^)</') is the radius of the NCS. We always set 
^min = -Rns for alTPC (and TPC) and R^,^ = i?Ncs 
for alOG (and alOG) fits, while R^^.^ may vary and 
may even be quite close to i?LC- Our alOG and alTPC 
models have 9 and 8 free parameters respectively, de- 
scribing the pulsar geometry (o; and C,) and gamma- 
ray and radio gap locations, apart from P which deter- 
mines the size of the PC (see Figure [1|). More details 
are provided in [511] . Note that the radio and gamma- 
ray emission layers are fit independently. 
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3.2. Low-altitude Gamma-ray and Radio 
Emission: laSG Models 

This model provides a non-caustic explanation for 
the emission, and may be viewed as a geometric low- 
altitude SG model [3^ resembling a hollow cone beam 
close to the stellar surface. We modulate the emissiv- 
ity according to 

fexpjA.Mn) .<.f 
[ exp(-As/crout) , s > Sf, 

with s the distance above the NS surface along a mag- 
netic field line, As = s — Sf , and din and CTout setting 
the rate at which the intensity rises and falls along 
the magnetic field lines. The peak intensity occurs at 
a distance s — (i.e.. As = 0) along the field lines. 
The laSG models have 5 free parameters (more details 

in mi). 



4. FINDING OPTIMAL LC FITS 

In order to statistically pick the best-fit parameters, 
for the alOG and alTPC models and for the three 
MSPs considered here, we have developed an Markov 
chain Monte Carlo (MCMC) maximum likelihood pro- 
cedure [S^l- The gamma-ray LCs are fit using Poisson 
likelihood while the radio LCs are fit using a statis- 
tic, and the two values are then combined. For a given 
parameter state the likelihood value is calculated by 
independently optimizing the radio and gamma-ray 
model normalizations using the scipy python mod- 
ule^ and the scipy. optimize. fmin_Lfbgs_b multivariate, 
bound optimizer [s^ . 

An MCMC involves taking random steps in param- 
eter space and accepting a step based on the likeli- 
hood ratio with respect to the previous step 27']. The 
likelihood surfaces can be very multimodal which can 
lead to poor mixing of the chain and slow conver- 
gence. Therefore, we have implemented small-world 
chain steps [l!] and simulated annealing [35.] to speed 
up the convergence and ensure that the MCMC fully 
explores the parameter space and does not get stuck 
in a local maximum. We verify that our chains have 
converged using the criteria proposed by [2l| . 

In order to balance the gamma-ray and radio con- 
tributions to the likelihood, we have chosen to use a 
radio uncertainty which is equal to the average rel- 
ative gamma-ray uncertainty in the on-peak region 
times the maximum radio value. The choice of ra- 
dio uncertainty can strongly affect the best-fit re- 
sults; in particular, a smaller uncertainty will decrease 
the overall likelihood and can, in some cases, lead to 



^ See |http: / /docs. scipy.org/ doc/ 1 for documentation. 
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Figure 2: LC fits for PSR J0034-0534 using alOG / 
alTPC models. Panel (a) shows the gamma-ray data, 
while panel (b) shows the radio data. 

a different best-fit geometry which favors the radio 
LC more strongly. When varying the radio uncer- 
tainty by a factor of 2, the best-fit a and ( values of 
PSR J0034-0534 were found to change by < 13°. For 
PSR J1939-h2134, the best-fit a and C were found to 
vary by < 7° when varying the radio uncertainty. The 
best-fit geometry of PSR J1959-h2048 was found to be 
the most sensitive to changes in the radio uncertainty, 
with either the best- fit a or ^ value changing by ~ 35°, 
while the other parameter changed by ^ 15°. 

Starting from the best-fit parameters found by the 
MCMC, we produced confidence contours in a and ^ 
by performing likelihood profile scans over the other 
parameters, allowing for the possibility of finding a 
better fit. The uncertainties on a and ^ quoted in Ta- 
bleware approximate 95% confidence level uncertain- 
ties. We can also estimate uncertainties on the emis- 
sion altitude parameters using the information from 
the likelihood profile scans which generated the confi- 
dence contours. Note that we used manually-selected 
LC fits for the laSG models. 



5. RESULTS 

The effects of letting i^min and i?max be free parame- 
ters in the alOG / alTPC model context, as well as us- 
ing different fading parameters in our laSG models is 
discussed elsewhere [5l|. Our best-fit LC parameters 
are summarized in TableU) As an example, the alOG / 
alTPC LC fits for PSR J0034-0534 are shown in Fig- 
ure [3 while Figure [3] shows fits for PSR J0034-0534 
in the case of laSG models. 



6. CONCLUSIONS 

We studied a third subclass of gamma-ray MSP LCs 
for which the gamma-ray and radio profiles are phase- 
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Figure 3: LC fits for PSR J0034-0534 using laSG 
models. Note that we had to introduce a large phase shift 
of 00 = 0.68 for both gamma-ray and radio LCs, as the 
data and model zero phases do not coincide. This implies 
that the definition of 'leading peak' and 'first peak' do 
not coincide anymore. See Table |I] for more details. 



aligned, and the gamma-ray and radio emission should 
therefore be co-located. We introduced free parame- 
ters i?min and i?max in the alOG / alTPC models, 
(both are free for the radio emission region, but only 
^max i^ fr^^ f*^^' the the gamma-ray emission region) 
and found fits from these models which could repro- 
duce the salient features of the profiles, although not 
perfectly. As a second option, we implemented the 
the laSG models and demonstrated that a modulated 
emissivity at low altitudes can reproduce main fea- 
tures of the profiles quite well. 

At the moment, it is difficult to quantitatively fa- 
vor one class of models above the other, since search- 
ing for optimal laSG LC fits has been done manu- 
ally. However, we did calculate the likelihood of the 
best-fit laSG LCs. The alTPC models provide slightly 
better LC fits than the alOG models, and both of 
these give better fits than the laSG models for the 
parameters listed (see Table U). Favoring the alOG / 
alTPC models over the laSG model therefore implies 
that the phase-aligned gamma-ray and radio LCs are 
most probably of caustic origin, produced in the outer 
magnetosphere, and the radio emission is most likely 
originating near the light cylinder. Thus, we can now 
divide the gamma-ray MSP population into three sub- 
classes on the basis of their LCs: those with LCs fit by 
standard OG / TPC models, those with phase-aligned 
LCs fit by alOG / alTPC or laSG models, and those 
with LCs fit by PSPC models. 

Radio polarization measurements can be used to 
give independent constraints on the pulsar viewing 
geometry [ssj, complementing the gamma-ray model 
fits, although the traditional RVM [40| is not valid for 
the alTPC or alOG models where the radio peaks are 
caustics. Furthermore, the RVM is not expected to 
yield good results in the case of radio cone beam emis- 



sion in MSPs, as these beams should suffer significant 
distortions due to retardation and aberration [9] . This 
may account for the generally poor or non-existent 
RVM fits of MSP polarization data. 

Caustic models predict rapid PA swings coupled 
with depolarization [iTj . since although the emission 
originates from a large range of altitudes and magnetic 
field orientations, it is restricted to a narrow phase 
interval to form the peaks. These features seem to 
be present in radio polarization measurements of the 
modeled MSPs [i^ 113, H^. Polarization signatures 
are important to help discriminate between models 
with caustic emission (such as occurs in alOG / alTPC 
models) and non-caustic emission (e.g., in the laSG 
model). 

Future studies include development of full radiation 
models which will be able to reproduce the multiwave- 
length LC shapes, polarization properties, as well as 
the energy-dependent behavior of the spectra of the 
gamma-ray MSPs. Ways to increase pair production 
also need to be found, which may include investigation 
of offset-PC dipole magnetic fields ^] and higher- 
multipole magnetic fields near the NS surface |5q . 



Table I Inferred best-fit model LC parameters for 

PSR J0034-0534, PSR J1939-H2134, and 

PSR J1959-f 2048. The columns represent the geometric 

model ('laSGl' refers to an laSG model with Sf = 1.27?, 

(Jin = Q.IR, and (Tout = 0.3-R, and 'laSG2' refers to an 

laSG model with Sf = 1.5ii, (Tin = 0.27i, and 

(Tout = 0.5/?), inclination and observer angles a and C, 

(measured in degrees), maximum gamma-ray altitude 

^max, minimum radio altitude -Rmim maximum radio 

altitude R^ax, as well as the log- likelihood 

A = — A In (like) of the fit. The altitudes are in units of 

Ri^c. We used R'^^^ = i?NCS for the alOG model, and 

■Rmiu = -Rns for the alTPC model. 



Model 


Q 


c 


m 

^ •'max 


-'•'min 




A 


PSR J0034 


0534 










alOG 
alTPC 


12+_f 
30l? 


mtf 

70±2 


0.9lg;? 
0.9±0.1 


2+° '^ 

"■^-0.06 
^■'-0.3 




96.1 
87.0 


laSGl 


10 


34 








97.3 


laSG2 


10 


37 








98.7 


PSR J1939+2134 


alOG 
alTPC 


84t^6 


84t^ 


1 Q+0.2 

1.0±0.2 


0.6±0.1 

' -0.3 


0.9±0.1 

0.91°:? 


130.9 
126.3 


laSGl 


30 


32 








146.9 


laSG2 


35 


25 








154.6 


PSR J1959+2048 


alOG 
alTPC 


31+f 
47t?3 


89l3 
85ti 


1 1+0.1 
^■^-0.2 

-'-•^-0.4 


0.7±0.1 
0.8±0.1 


0.9+°:? 
i-ot^:? 


128.3 
123.7 


laSGl 


20 


43 








129.0 


laSG2 


25 


45 








141.7 
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